function [simul,pr_entrant] = DO_sim_func(simul_size,param,VO,VN,U,U_hat,opt_ind_hp_N,opt_ind_hp_O,GU,GN,GO,opt_p_U,opt_p_N,opt_p_O,opt_x_N,opt_x_O,opt_eta_U,opt_eta_N,opt_eta_O)

    %----- Specify Parameters -----%
    nu = param(1);
    beta = param(2);            % Discount factor (factors in exit probability)

    d = param(4);               % Human capital depreciation rate
    lambda = param(6);          % Probability of search OTJ
    chi = param(7);             % Investment cost parameter

    delta =param(11);           % Separation rate
    h_sigma =param(12);         % Standard deviation of initial human capital distribution

    pr_phi(1) = param(13);      % Probability of entering the model with phi=1 (low fixed productivity type)
    pr_phi(2) = param(14);      % Probability of entering the model with phi=2 (high fixed productivity type)
    gamma=2;

    % Specify grids %
    hlb = 0.0;                 % Human capital grid lower bound 
    hub = 7.0;                 % Human capital grid upper bound 
    nh = 3000;                 % Number of possible human capital values
    h_grid = linspace(hlb,hub,nh); % Create human capital grid
    
    % ======================= Create discretized h distribution ==================== %
    h_mu = 1;                  % Mean of initial human capital distribution
    cdf_h = zeros(1,nh);       % CDF of initial h distribution
    pdf_h = zeros(1,nh);       % PDF of initial h distribution

    h_bin_size = (hub-hlb)/(nh-1);
    for ind_h = 1:nh
        if (ind_h ==1)
            bin_ub_now = h_grid(ind_h) + (h_bin_size)/2;
            cdf_h(ind_h) = normcdf(bin_ub_now,h_mu,h_sigma);
            pdf_h(ind_h) = normcdf(bin_ub_now,h_mu,h_sigma);
        else
            bin_ub_now = h_grid(ind_h) + (h_bin_size)/2;
            bin_ub_last = h_grid(ind_h-1) + (h_bin_size)/2;
            cdf_h(ind_h) = normcdf(bin_ub_now,h_mu,h_sigma);
            pdf_h(ind_h) = normcdf(bin_ub_now,h_mu,h_sigma) - normcdf(bin_ub_last,h_mu,h_sigma);
        end
    end

    % Find index for next period human capital (h') if the agent only lets human capital depreciate
    ind_hp_min = ones(1,nh);          % Minimum index for h' (index for h' if human capital only depreciates)
    for ind_h = 1:nh
        h = h_grid(ind_h);
        error_min = 100;
        for ind_hp = 1:nh
            hp = h_grid(ind_hp);       % hp is the value of h'
            error = abs(hp - h*(1-d)); % Find index that minimizes (hp-h*(1-d))
            if (error<error_min)
                error_min = error;
                ind_hp_min(ind_h) = ind_hp;
            end
        end
    end
    
    %% =============================================================================================================== %%
    %======================================== Simulation for Wage Moments and Results ==================================%
    %================================================================================================================= %%
    
    % Need to follow "respondents" from entrance because SS distribution does not save the eta values %
    nt = 400;               % Number of quarters to follow each individual - 100 yrs max - most will exit the model before this
    simul(simul_size,nt) = struct;
    rng(123456789);         % Set random number seed so results are exactly the same each time the code is ran   

    % Aso want to keep track of transition probabilities for entrants
    sum_entrant_O = 0;
    sum_entrant_N = 0;
    sum_entrant_U = 0;
    
    for id = 1:simul_size   % Loop over each individual in the sample
       
        for t = 1:nt        % Loop over each time period the individual is observed
            
            if (t==1)
                % For new entrants, assign phi and h, then see if they are successful in their job search from unemployment %
                % Assign starting h value
                r1 = rand;
                simul(id,t).h_ind = 1;  % If r1 is not bigger than any CDF values, assign to the lowest slot
                for ind_h = 2:nh
                    if (r1>cdf_h(ind_h-1)) 
                        simul(id,t).h_ind = ind_h;
                    end 
                end
                ind_h = simul(id,t).h_ind;
                h = h_grid(ind_h);
                % Now assign phi type
                r1 = rand;
                if r1 < pr_phi(1)
                    simul(id,t).phi_ind = 1;
                else
                    simul(id,t).phi_ind = 2;
                end
                ind_phi = simul(id,t).phi_ind;
                % See if new entrant is successful in job search from unemployment - - If successful, save their wage %
                r1 = rand;
                if (r1 < opt_p_U(ind_phi,ind_h)) % Finds job
                    if (GU(ind_phi,ind_h)==1 || GU(ind_phi,ind_h)==3)        % Search for N
                        simul(id,t).stat = 1;                                % Status "stat" 1 = in non-outsourced job
                        sum_entrant_N = sum_entrant_N + 1;
                        simul(id,t).eta = opt_eta_U(ind_phi,ind_h);
                        eta = simul(id,t).eta;
                        
                        ind_hp = opt_ind_hp_N(ind_phi,ind_h);                % This is the h' value they will choose for next period 
                        hp = h_grid(ind_hp);
                        simul(id,t).wage = eta*(VN(ind_phi,ind_h) - U(ind_phi,ind_h)) + U(ind_phi,ind_h) + (chi*(hp - h*(1-d)))^gamma...
                            - beta*delta*U_hat(ind_phi,ind_hp) - beta*(1-delta)*lambda*opt_p_N(ind_phi,ind_hp)*opt_x_N(ind_phi,ind_hp)...
                            - beta*(1-delta)*(1-lambda*opt_p_N(ind_phi,ind_hp))*(eta*VN(ind_phi,ind_hp) + (1-eta)*U(ind_phi,ind_hp));
                            
                    else                                                      % Search for O
                        simul(id,t).stat = 2;                                 % Status "stat" 2 = in outsourced job 
                        sum_entrant_O = sum_entrant_O + 1;
                        simul(id,t).eta = opt_eta_U(ind_phi,ind_h);
                        eta = simul(id,t).eta;
                        
                        ind_hp = opt_ind_hp_O(ind_phi,ind_h);                % This is the h' value they will choose for next period 
                        hp = h_grid(ind_hp);
     
                        simul(id,t).wage = eta*(VO(ind_phi,ind_h) - U(ind_phi,ind_h)) + U(ind_phi,ind_h) + (chi*(hp - h*(1-d)))^gamma...
                            - beta*delta*U_hat(ind_phi,ind_hp) - beta*(1-delta)*lambda*opt_p_O(ind_phi,ind_hp)*opt_x_O(ind_phi,ind_hp)...
                            - beta*(1-delta)*(1-lambda*opt_p_O(ind_phi,ind_hp))*(eta*VO(ind_phi,ind_hp) + (1-eta)*U(ind_phi,ind_hp));   
                    end   
                else  % Does not find job
                    simul(id,t).stat = 0;                                    % Status "stat" 0 = unemployed  
                    sum_entrant_U = sum_entrant_U + 1;
                    simul(id,t).wage = 0.0;
                    simul(id,t).eta = 0.0;
                end   % End if t=1 individual finds job or not
                
            else % so t>1
                simul(id,t).phi_ind = simul(id,t-1).phi_ind;                 % Individual's phi value is permanent
                ind_phi = simul(id,t).phi_ind;
                
                if (simul(id,t-1).stat==999) % if dead                       % If the agent already exited the model they remain dead             
                    simul(id,t).stat=999;                                    % Status "stat" 999 = dead
                else   % was not dead in last period
                    rnu = rand;
                    if (rnu<nu) % Dies at end of last period
                        simul(id,t).stat=999;
                    else        % Does NOT die at end of last period
                        
                        if ((simul(id,t-1).stat)==0)        % If unemployed last period 
                            % Recall h from last period and apply any depreciation
                            ind_h = simul(id,t-1).h_ind;
                            simul(id,t).h_ind =  ind_hp_min(ind_h);
                            ind_h = simul(id,t).h_ind;
                            h = h_grid(ind_h);
                            
                            % If unemployed last period - see if matches from unemployment !
                            r1 = rand;
                            if (r1 < opt_p_U(ind_phi,ind_h))    % Finds job
                                if (GU(ind_phi,ind_h)==1 || GU(ind_phi,ind_h)==3)      % Search for N job
                                    simul(id,t).stat = 1; 
                                    simul(id,t).eta = opt_eta_U(ind_phi,ind_h);
                                    eta = simul(id,t).eta;
                                    ind_hp = opt_ind_hp_N(ind_phi,ind_h);
                                    hp = h_grid(ind_hp);
                                    simul(id,t).wage = eta*(VN(ind_phi,ind_h) - U(ind_phi,ind_h)) + U(ind_phi,ind_h) + (chi*(hp - h*(1-d)))^gamma...
                                        - beta*delta*U_hat(ind_phi,ind_hp) - beta*(1-delta)*lambda*opt_p_N(ind_phi,ind_hp)*opt_x_N(ind_phi,ind_hp)...
                                        - beta*(1-delta)*(1-lambda*opt_p_N(ind_phi,ind_hp))*(eta*VN(ind_phi,ind_hp) + (1-eta)*U(ind_phi,ind_hp));  
                                else                                                    % Search for O job
                                    simul(id,t).stat = 2; 
                                    simul(id,t).eta = opt_eta_U(ind_phi,ind_h);
                                    eta = simul(id,t).eta;
                                    ind_hp = opt_ind_hp_O(ind_phi,ind_h);
                                    hp = h_grid(ind_hp);
                                    simul(id,t).wage = eta*(VO(ind_phi,ind_h) - U(ind_phi,ind_h)) + U(ind_phi,ind_h) + (chi*(hp - h*(1-d)))^gamma...
                                        - beta*delta*U_hat(ind_phi,ind_hp) - beta*(1-delta)*lambda*opt_p_O(ind_phi,ind_hp)*opt_x_O(ind_phi,ind_hp)...
                                        - beta*(1-delta)*(1-lambda*opt_p_O(ind_phi,ind_hp))*(eta*VO(ind_phi,ind_hp) + (1-eta)*U(ind_phi,ind_hp));  
                                end
                            else                                 % Unemployed last period and does not find job%
                                simul(id,t).stat = 0;
                                simul(id,t).wage = 0.0;
                                simul(id,t).eta = 0.0;
                            end
     
                        elseif ((simul(id,t-1).stat)==1)    % If in N last period  
                            % Recall h from last period and apply any investment
                            ind_h = simul(id,t-1).h_ind;
                            simul(id,t).h_ind = opt_ind_hp_N(ind_phi,ind_h);
                            ind_h = simul(id,t).h_ind;
                            h = h_grid(ind_h);
                            
                            % If N last period - see if separation shock %
                            rdelta = rand;
                            if (rdelta < delta)  % N last period - separation shock    %
                                % See if finds job from unemployment
                                r1 = rand;
                                if (r1 < opt_p_U(ind_phi,ind_h)) % N last period - separation shock - finds job from unemployment
                                    if (GU(ind_phi,ind_h)==1 || GU(ind_phi,ind_h)==3)      % Search for N
                                        simul(id,t).stat = 1; 
                                        simul(id,t).eta = opt_eta_U(ind_phi,ind_h);
                                        eta = simul(id,t).eta;
                                        ind_hp = opt_ind_hp_N(ind_phi,ind_h);
                                        hp = h_grid(ind_hp);
                                        simul(id,t).wage = eta*(VN(ind_phi,ind_h) - U(ind_phi,ind_h)) + U(ind_phi,ind_h) + (chi*(hp - h*(1-d)))^gamma...
                                            - beta*delta*U_hat(ind_phi,ind_hp) - beta*(1-delta)*lambda*opt_p_N(ind_phi,ind_hp)*opt_x_N(ind_phi,ind_hp)...
                                            - beta*(1-delta)*(1-lambda*opt_p_N(ind_phi,ind_hp))*(eta*VN(ind_phi,ind_hp) + (1-eta)*U(ind_phi,ind_hp));  
                                    else                                                    % Search for O
                                        simul(id,t).stat = 2; 
                                        simul(id,t).eta = opt_eta_U(ind_phi,ind_h);
                                        eta = simul(id,t).eta;
                                        ind_hp = opt_ind_hp_O(ind_phi,ind_h);
                                        hp = h_grid(ind_hp);
                                        simul(id,t).wage = eta*(VO(ind_phi,ind_h) - U(ind_phi,ind_h)) + U(ind_phi,ind_h) + (chi*(hp - h*(1-d)))^gamma...
                                            - beta*delta*U_hat(ind_phi,ind_hp) - beta*(1-delta)*lambda*opt_p_O(ind_phi,ind_hp)*opt_x_O(ind_phi,ind_hp)...
                                            - beta*(1-delta)*(1-lambda*opt_p_O(ind_phi,ind_hp))*(eta*VO(ind_phi,ind_hp) + (1-eta)*U(ind_phi,ind_hp)); 
                                    end
                                else                                 % N last period - separation shock - does find NOT job from unemployment
                                    simul(id,t).stat = 0;
                                    simul(id,t).wage = 0.0;
                                    simul(id,t).eta = 0.0;
                                end
                        
                            else                   % N last period - no separation shock %
                                % See if successfully searchs OTJ
                                r1 = rand;
                                if (r1 < lambda*opt_p_N(ind_phi,ind_h)) % N last period - no separation shock - successfully searchs OTJ
                            
                                    if (GN(ind_phi,ind_h)==1 || GN(ind_phi,ind_h)==3)             % N - Search OTJ for N
                                        simul(id,t).stat = 1; 
                                        simul(id,t).eta = opt_eta_N(ind_phi,ind_h);
                                        eta = simul(id,t).eta;
                                        ind_hp = opt_ind_hp_N(ind_phi,ind_h);
                                        hp = h_grid(ind_hp);
                                        simul(id,t).wage = eta*(VN(ind_phi,ind_h) - U(ind_phi,ind_h)) + U(ind_phi,ind_h) + (chi*(hp - h*(1-d)))^gamma...
                                            - beta*delta*U_hat(ind_phi,ind_hp) - beta*(1-delta)*lambda*opt_p_N(ind_phi,ind_hp)*opt_x_N(ind_phi,ind_hp)...
                                            - beta*(1-delta)*(1-lambda*opt_p_N(ind_phi,ind_hp))*(eta*VN(ind_phi,ind_hp) + (1-eta)*U(ind_phi,ind_hp)); 
                                    else                                                           % N - Search OTJ for O
                                        simul(id,t).stat = 2; 
                                        simul(id,t).eta = opt_eta_N(ind_phi,ind_h);
                                        eta = simul(id,t).eta;
                                        ind_hp = opt_ind_hp_O(ind_phi,ind_h);
                                        hp = h_grid(ind_hp);
                                        simul(id,t).wage = eta*(VO(ind_phi,ind_h) - U(ind_phi,ind_h)) + U(ind_phi,ind_h) + (chi*(hp - h*(1-d)))^gamma...
                                            - beta*delta*U_hat(ind_phi,ind_hp) - beta*(1-delta)*lambda*opt_p_O(ind_phi,ind_hp)*opt_x_O(ind_phi,ind_hp)...
                                            - beta*(1-delta)*(1-lambda*opt_p_O(ind_phi,ind_hp))*(eta*VO(ind_phi,ind_hp) + (1-eta)*U(ind_phi,ind_hp));    
                                    end
                                else                                         % N last period - no separation shock - does NOT successfully search OTJ 
                                    simul(id,t).stat = simul(id,t-1).stat;
                                    simul(id,t).eta = simul(id,t-1).eta;
                                    eta = simul(id,t).eta;
                                    ind_hp = opt_ind_hp_N(ind_phi,ind_h);
                                    hp = h_grid(ind_hp);
                                    simul(id,t).wage = eta*(VN(ind_phi,ind_h) - U(ind_phi,ind_h)) + U(ind_phi,ind_h) + (chi*(hp - h*(1-d)))^gamma...
                                        - beta*delta*U_hat(ind_phi,ind_hp) - beta*(1-delta)*lambda*opt_p_N(ind_phi,ind_hp)*opt_x_N(ind_phi,ind_hp)...
                                        - beta*(1-delta)*(1-lambda*opt_p_N(ind_phi,ind_hp))*(eta*VN(ind_phi,ind_hp) + (1-eta)*U(ind_phi,ind_hp));
                                end
                            end % End if receives separation shock when N last period     
                        else                                                  % If in O last period  
                            % Recall h from last period and apply any investment
                            ind_h = simul(id,t-1).h_ind; 
                            simul(id,t).h_ind = opt_ind_hp_O(ind_phi,ind_h);
                            ind_h = simul(id,t).h_ind;
                            h = h_grid(ind_h);
                            
                            % If O last period - see if separation shock %
                            rdelta = rand;
                            if (rdelta < delta) % O last period - separation shock    
                                % See if finds job from unemployment
                                r1 = rand;
                                if (r1 < opt_p_U(ind_phi,ind_h))    % N last period - separation shock - finds job from unemployment
                                    if (GU(ind_phi,ind_h)==1 || GU(ind_phi,ind_h)==3)       % Search for N
                                        simul(id,t).stat = 1; 
                                        simul(id,t).eta = opt_eta_U(ind_phi,ind_h);
                                        eta = simul(id,t).eta;
                                        ind_hp = opt_ind_hp_N(ind_phi,ind_h);
                                        hp = h_grid(ind_hp);
                                        simul(id,t).wage = eta*(VN(ind_phi,ind_h) - U(ind_phi,ind_h)) + U(ind_phi,ind_h) + (chi*(hp - h*(1-d)))^gamma...
                                            - beta*delta*U_hat(ind_phi,ind_hp) - beta*(1-delta)*lambda*opt_p_N(ind_phi,ind_hp)*opt_x_N(ind_phi,ind_hp)...
                                            - beta*(1-delta)*(1-lambda*opt_p_N(ind_phi,ind_hp))*(eta*VN(ind_phi,ind_hp) + (1-eta)*U(ind_phi,ind_hp));
                                    else                                                     % Search for O
                                        simul(id,t).stat = 2; 
                                        simul(id,t).eta = opt_eta_U(ind_phi,ind_h);
                                        eta = simul(id,t).eta;
                                        ind_hp = opt_ind_hp_O(ind_phi,ind_h);
                                        hp = h_grid(ind_hp);
                                        simul(id,t).wage = eta*(VO(ind_phi,ind_h) - U(ind_phi,ind_h)) + U(ind_phi,ind_h) + (chi*(hp - h*(1-d)))^gamma...
                                            - beta*delta*U_hat(ind_phi,ind_hp) - beta*(1-delta)*lambda*opt_p_O(ind_phi,ind_hp)*opt_x_O(ind_phi,ind_hp)...
                                            - beta*(1-delta)*(1-lambda*opt_p_O(ind_phi,ind_hp))*(eta*VO(ind_phi,ind_hp) + (1-eta)*U(ind_phi,ind_hp)); 
                                    end
                                
                                else                                 % O last period - separation shock - does find NOT job from unemployment
                                    simul(id,t).stat = 0;
                                    simul(id,t).wage = 0.0;
                                    simul(id,t).eta = 0.0;
                                end
                            else                   % O last period - no separation shock %
                                %See if successfully searchs OTJ
                                r1 = rand;
                                if (r1 < lambda*opt_p_O(ind_phi,ind_h))   % O last period - no separation shock - successfully searchs OTJ
                                    if (GO(ind_phi,ind_h)==1 || GO(ind_phi,ind_h)==3)            % O - Search OTJ for N
                                        simul(id,t).stat = 1; 
                                        simul(id,t).eta = opt_eta_O(ind_phi,ind_h);
                                        eta = simul(id,t).eta;
                                        ind_hp = opt_ind_hp_N(ind_phi,ind_h);
                                        hp = h_grid(ind_hp);
                                        simul(id,t).wage = eta*(VN(ind_phi,ind_h) - U(ind_phi,ind_h)) + U(ind_phi,ind_h) + (chi*(hp - h*(1-d)))^gamma...
                                            - beta*delta*U_hat(ind_phi,ind_hp) - beta*(1-delta)*lambda*opt_p_N(ind_phi,ind_hp)*opt_x_N(ind_phi,ind_hp)...
                                            - beta*(1-delta)*(1-lambda*opt_p_N(ind_phi,ind_hp))*(eta*VN(ind_phi,ind_hp) + (1-eta)*U(ind_phi,ind_hp));  
                                    else                                                          % O - Search OTJ for O
                                        simul(id,t).stat = 2; 
                                        simul(id,t).eta = opt_eta_O(ind_phi,ind_h);
                                        eta = simul(id,t).eta;
                                        ind_hp = opt_ind_hp_O(ind_phi,ind_h);
                                        hp = h_grid(ind_hp);
                                        simul(id,t).wage = eta*(VO(ind_phi,ind_h) - U(ind_phi,ind_h)) + U(ind_phi,ind_h) + (chi*(hp - h*(1-d)))^gamma...
                                            - beta*delta*U_hat(ind_phi,ind_hp) - beta*(1-delta)*lambda*opt_p_O(ind_phi,ind_hp)*opt_x_O(ind_phi,ind_hp)...
                                            - beta*(1-delta)*(1-lambda*opt_p_O(ind_phi,ind_hp))*(eta*VO(ind_phi,ind_hp) + (1-eta)*U(ind_phi,ind_hp));  
                                    end
                                else                                         % O last period - no separation shock - does NOT successfully searchs OTJ
                                    simul(id,t).stat = simul(id,t-1).stat;
                                    simul(id,t).eta = simul(id,t-1).eta;
                                    eta = simul(id,t).eta;
                                    ind_hp = opt_ind_hp_O(ind_phi,ind_h);  
                                    hp = h_grid(ind_hp);
                                    simul(id,t).wage = eta*(VO(ind_phi,ind_h) - U(ind_phi,ind_h)) + U(ind_phi,ind_h) + (chi*(hp - h*(1-d)))^gamma...
                                        - beta*delta*U_hat(ind_phi,ind_hp) - beta*(1-delta)*lambda*opt_p_O(ind_phi,ind_hp)*opt_x_O(ind_phi,ind_hp)...
                                        - beta*(1-delta)*(1-lambda*opt_p_O(ind_phi,ind_hp))*(eta*VO(ind_phi,ind_hp) + (1-eta)*U(ind_phi,ind_hp));  
                                end
                            end % End if receives separation shock when OL last period
                        end % End if for stat (when not dead or died before period start) last period
                    end         % End if dies at end of last period
                end    % End if dead in last period
            end % End if t==1
        end % End loop over time periods in panel
    end % End loop over ids in simul_size
    
    pr_entrant_N = sum_entrant_N/simul_size;
    pr_entrant_O = sum_entrant_O/simul_size;
    pr_entrant_U = sum_entrant_U/simul_size;

    pr_entrant = [pr_entrant_U,pr_entrant_N,pr_entrant_O];

end % End function